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Analysis of linearized Galerkin-mixed FEMs for the 
time-dependent Ginzbnrg-Landau eqnations of snpercondnctivity 


Huadong Gao * and Weiwei Sun ^ 


Abstract 

A linearized backward Euler Galerkin-mixed finite element method is investigated for 
the time-dependent Ginzbnrg-Landau (TDGL) equations under the Lorentz gauge. By in¬ 
troducing the induced magnetic field cr = curl A as a new variable, the Galerkin-mixed 
FE scheme offers many advantages over conventional Lagrange type Galer kin FEMs. An 
optimal error estimate for the linearized Galerkin-mixed FE scheme is established uncondi¬ 
tionally. Analysis is given under more general assumptions for the regularity of the solution 
of the TDGL equations, which includes the problem in two-dimensional noncovex polygons 
and certain three dimensional polyhedrons, while the conventional Galerkin FEMs may not 
converge to a true solution in these cases. Numerical examples in both two and three di¬ 
mensional spaces are presented to confirm our theoretical analysis. Numerical results show 
clearly the efficiency of the mixed method, particularly for problems on nonconvex domains. 

Keywords: Ginzburg-Landau equation, linearized scheme, mixed finite element method, 
unconditional convergence, optimal error estimate, superconductivity. 


1 Introduction 

In this paper, we consider the time-dependent Ginzburg-Landau (TDGL) equations under the 
Lorentz gauge 

dll) % 

rj^ - ir]K{div A)ijj + (-V -b A)^V’ + (IV'I^ - 1)V’ = 0 , (1.1) 

UL K 

^/V 9 

—-Vdiv A + curl curl A H- — ipV'ijj*) -|- IV'I^A = curl Hg , (1.2) 

(yt> 2 av 

for X E and t E (0, T], where the complex scalar function if) is the order parameter and the real 
vector-valued function A is the magnetic potential. In (ll.ip - (ll.2D . denotes the density of the 
superconducting electron pairs. |pp = 1 and = 0 represent the perfectly superconducting 
state and the normal state, respectively, while 0 < < 1 represents a mixed state. The 

real vector-valued function Hg is the external applied magnetic field, k is the Ginzburg-Landau 
parameter and rj is a dimensionless constant. In the following, we set rj = 1 for the sake of 
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simplicity. We assume that 0 is a simply-connected bounded Lipschitz domain in The 
following boundary and initial conditions are supplemented to ([nD-([L2D 

dib 

= 0, curl A X n = He x n, A • n = 0, on (90 x [0,T], (1-3) 

A(x,0) = Ao(x), inO, (1.4) 

where n is the outward unit normal vector. We note that in (ll.ip - (ll.4p . V, div and curl are the 
standard gradient, divergence and rotation operators in three dimensional space. 

The TDGL equations were first deduced by Gor’kov and Eliashberg in [26] from the mi¬ 
croscopic Bardeen-Cooper-Schrieffer theory of superconductivity. For detailed physical de¬ 
scription and mathematical modeling of the superconductivity phenomena, we refer to the re¬ 
view articles uniEQ]. Theoretical analysis for the TDGL equations can be found in literature 
[IllSIllMllMlEHl. The global existence and uniqueness of the strong solution were established in 
[14| for the TDGL equations with the Lorentz gauge. Numerical methods for solving the TDGL 
equations have also been studied extensively, e.g., see lainillsllIllllsllSTllMllMlEilSTlllI]. A 
semi-implicit Euler scheme with a finite element approximation was first proposed by Chen and 
Hoffmann m for the TDGL equations with the Lorentz gauge. A suboptimal error estimate 
was obtained for the equations in two dimensional space. Later, a decoupled alternating Crank- 
Nicolson Galerkin method was proposed by Mu and Huang |35j . An optimal error estimate was 
presented under the time step restrictive conditions r = 0(/i 12 ) for the two-dimensional model 
and T = 0{h?‘) for the three-dimensional model, where h and r are the mesh size in the spatial 
direction and the time direction, respectively. All these schemes in HaESllMI are nonlinear. At 
each time step, one has to solve a nonlinear system. Clearly, more commonly-used discretiza¬ 
tions for nonlinear parabolic equations are linearized schemes, which only require the solution 
of a linear system at each time step. A linearized Crank-Nicolson type scheme was proposed in 
|34j for slightly different TDGL equations and a systematic numerical investigation was made 
there. An optimal error estimate of a linearized Crank-Nicolson Galerkin EE scheme for the 
TDGL equations was provided unconditionally in a recent work |23j . 

All these methods mentioned above were introduced to solve the TDGL equations (|1.1I1 - (I1.4I) 
for the order parameter ijj and the magnetic potential A by Galerkin FEMs (or finite difference 
methods) and then, to calculate the magnetic field curl A by certain numerical differentiation. 
There are several drawbacks for such approaches. One of important issues in the study of the 
vortex motion in superconductors is the influence of geometric defects, which is of high interest in 
physics and can be considered as a problem in polygons. Previous numerical results [22l|23l|34| 
by conventional Lagrange FEMs showed the singularity and inaccuracy of the numerical solution 
around corners of polygons. Moreover, for the problem in a nonconvex polygon, conventional 
Lagrange FEMs for the TDGL equations may converge to a spurious solution, see the numerical 
experiments reported in |22l [32] for the problem on an L-shape domain. Another issue in the 
superconductivity model is its coupled boundary conditions in (II.3h . which may bring some 
extra difficulties in implementation to conventional Lagrange FEMs on a general domain, see 
section 5 of m for more detailed description. It is natural that a mixed finite element method 
with the physical quantity cr = curl A being the new variable may offer many advantages. 

It has been noted that mixed methods for second order elliptic problems, such as the scalar 
Poisson equation —Art = / with simple boundary conditions, Dirichlet or Neumann, have been 
well studied with the extra variable u = Vu in the last several decades, e.g., see PEllMllSSlIio] 
and references therein. However, mixed finite element methods for the vector Poisson equation 
— Au = f with the coupled boundary conditions curlu x n = g and u • n = 0 seems much 
more complicated than for a scalar equation. Theoretical analysis was established quite recently 
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in terms of finite element exterior calculus. For the TDGL equations (I2.13p - (I2.16I) in 
two dimensional space, Chen m proposed a semi-implicit and weakly nonlinear scheme with 
a mixed finite element method, in which both curl A and divA were introduced to be extra 
unknowns, and the corresponding finite element spaces for curl A and divA were constructed 
with certain bubble functions and some special treatments in elements adjacent to the boundary 
of the domain were also needed. A suboptimal error estimate was provided and no numerical 
example was given in [12] . Recently, a linearized backward Euler Galerkin-mixed finite element 
method was proposed in [22] for the TDGL equations in both two and three dimensional spaces 
with only one extra physical quantity a = curl A (<7 = curl A in two dimensional space). No 
analysis was provided in [22|. The scheme is decoupled and at each time step, one only needs 
to solve two linear systems for ijj and (A, cr), which can be done simultaneously. 

This paper focuses on theoretical analysis on a linearized Galerkin-mixed FEM for the TDGL 
equations to establish optimal error estimates. More important is that our analysis is given 
without any time-step restrictions and under more general assumptions for the regularity of 
the solution of the TDGL equations, which includes the problem in nonconvex polygons and 
certain convex polyhedrons. Usually, conventional EEMs for a scalar parabolic equation require 
the regularity of the solution in with s > 0 mm. while for the TDGL equations in 

a nonconvex polygon, A G H® and curl A, divA G with s < 1 in general [H]. Our 

numerical results show clearly that the mixed method converges to a true solution for problems 
in a nonconvex polygon and conventional Lagrange type FEMs do not in this case. In addition, 
in the mixed EEM, the coupled boundary condition curl A x n = He X n reduces to a Dirichlet 
type boundary condition <t x n = He x n. Implementation of Dirichlet boundary conditions for 
vector-valued elements (such as Nedelec and Raviart-Thomas elements) becomes much simpler. 
In a recent work by Li and Zhang, an approach based on Hodge decomposition was proposed 
and analyzed with optimal error estimates. However, the approach is applicable only for the 
problem in two-dimensional space due to the restriction of the Hodge decomposition. 

The rest of this paper is organized as follows. In section |2l we introduce the linearized back¬ 
ward Euler scheme with Galerkin-mixed finite element approximations for the TDGL equations 
and we present our main results on unconditionally optimal error estimates. In section jS] we 
recall some results for an auxiliary elliptic problem of the vector Poisson equation with the 
coupled boundary conditions. In section (H we prove that an optimal error estimate holds 
almost unconditionally {i.e., without any mesh ratio restriction). In section[5l we provide several 
numerical examples to confirm our theoretical analysis and show the efficiency of the proposed 
methods. Some concluding remarks are given in section [6l 

2 A linearized backward Euler Galerkin-mixed FEM 

In this section, we present a linearized backward Euler Galerkin-mixed finite element method for 
the TDGL equations and our main results. For simplicity, we introduce some standard notations 
and operators below. For any two complex functions u, v £ £^(D), we denote the T^(D) inner 
product and norm by 


{u,v) = / u{x) {v{x))* dx, \\u\\i ^2 = {u,u)2 , 

Jn 

where v* denotes the conjugate of the complex function v. Let IU^’P(D) be the Sobolev space 
defined on D, and by conventional notations, H^(D) := lU^’^(D), H^{Q) := lU^’^(D). Let 
'R^(D) = {u -|- iv\u,v G iL^(D)} be a complex-valued Sobolev space and H^(D) = [i/^(D)]'^ 
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be a vector-valued Sobolev space, where d is the dimension of fl. For a positive real number 
s = k + 9 with 0 < 0 < 1, we define H^{Q) = by the complex interpolation, see 

[6]. To introduce the mixed variational formulation, we denote 


H(div) = {A I A G L^(n), div A G with ||A||H(div) = (l|A |||2 + ||div AH^a) ^ 

and 

H(curl) = {a I A G L^(n), curl A G L^(0)} with || A||H(curi) = (II A |||2 + ||curl AIII 2 ) ^ • 
Also we dehne 


H(div) = {A I A G H(div), A • n 

o 

and its dual space H(div)' with norm 

(v, w) 


an 


l^llH(div)-= 

wGH(div) 


W 


H(div) 


0} 


Moreover, we denote 


H(curl) = IA I A G H(curl), A x = O} . 

By introducing cr = curl A, the mixed form of the TDGL equations (in])-([IlD can be 
written by 


^ - m(div A)'0 -I- (-V -|- A)‘^ip + (IV’P - 1)^’ = 0 , 
at K 

cr = curlA, 

^ yV a 

— -Vdiv A -I- curler -|-('0*VV’ — ip'S/'ijj*) + IV’pA = curl He , 

ut ‘2 iK 


with boundary and initial conditions 

dip 


dn 


= 0, <T X n = He X n, A • n = 0, 


on dQ X [0, T], 


V'(x,0) = ipo{x), cr{x,0) = curlAo(x), A(x,0) = Ao(x), in Q. 


( 2 . 1 ) 

( 2 . 2 ) 

(2.3) 


(2.4) 

(2.5) 


The mixed variational formulation of the TDGL equations ()2.ip - (l2.3l) with boundary and initial 
conditions (12.4^ - 112.5p is to find ip G L^(0, T;'H^(D)) with ^ G L^(0,T;'H“^(D)), and (cr, A) G 
L^(0, T; H(curl)) X L^(0, T; H(div)) with ^ G L^(0, T; H(div)'), where cr x n = He x n on dO., 
such that 

(^,w) - iK{{d\v A)ip,uj) + ((-V -h A)ip, {-V + A)a;) 

(JL Kj Kj 

+ {{\ip\^-i)ip,uj) = 0, yuj^n^n) (2.6) 

and 

(cr,x) - (curlx. A) = 0, V x G H(curl), (2.7) 

(^^, v) (curler, v) (div A,divv) -h —{i'ip*'^'ip - ipVip*),^) 
ut 2f^ 

-I-(IV^pA, v) = (curlHe, v), VvGH(div), (2.8) 
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for a.e. t G (0,T] with il!{x,0) = i/joix), cr{x,0) = curl Ao(a;) and A(x,0) = Ao(x). 

For simplicity, we assume that 0 is a polyhedron in three dimensional space. Let 7^ be a 
quasi-uniform tetrahedral partition of n with = Dk^k and denote hy h = max^j^e 7 ^{diam Qk} 
the mesh size. For a given partition 7h, we denote by the r-th order Lagrange finite element 
subspace of we denote by the r-th order first type Nedelec finite element subspace 

of H(curl), where the case r = 1 corresponds to the lowest order Nedelec edge element (6 dofs). 
We denote by the r-th order Raviart-Thomas finite element subspace of H(div), where the 
case r = 0 corresponds to the lowest order Raviart-Thomas face element (4 dofs). We also define 

o o o o 

Q)) = Q))nH(curl) and U)) = U))nH(div). It should be remarked that Brezzi-Douglas-Marini 
element can also be used for approximation of H(div). Here we only confine our attention to the 
Raviart-Thomas element approximation. By noting the approximation properties of the finite 
element spaces Q)) and U)) [HEllJl], we denote by vr^ a general projection operator on 
Qj) and Uj), satisfying 

||a; — TThUj\\j ^2 < Ch^\\co\\fjs , 0 < s < r -|- 1, 

< \\x-'^hX\\L^ + \\curl{x-7rhX)\\L^ <Ch^{\\x\\Ho + \\^^^^x\\Ho)^ l<s<r, (2.9) 

, ||r; - vr/iu ||^2 -h ||div {v - 7rhv)\\^2 < Ch^{\\v\\jjs + ||divu||jys), 0 < s < r -h 1. 

Let {tn}n=o ^ uniform partition in the time direction with the step size t = and let 
vP = u{-,nT). For a sequence of functions {C/"'}))Lq defined on Q, we denote 

jjn _ jjn-l 

DrU^ = -, forn = 1, 2, ..., A. 

r 

With the above notations, the linearized backward Euler Galerkin-mixed FEM for the mixed 
form TDGL equations (|2.1l) - (|2.5p is to find G and (<t)) , A))) G x Uj), with x n = 
X n on dQ, such that for n = 1,2,... ,N, 

{D,rh,o^u) - K~")rh,o^h) + ((-V + Al-^Wh , (-V + Al-^)u:n) 

tx K 

- l)i^^,uJh) = 0, Vo;;, G vL (2.10) 

and 

{(Xh,Xh) - icurlxh,-^l) = 0, VxfeGQ))+S (2.11) 

{DrAl, Vh) + (div Al, div Vh) + (curl , Vh) + (| A)(, v/j) 

= (curlH^v,)-^((V^^^)*V^^l-V’r'V(V’r')^v/.) , Vv.GU);, (2.12) 

where r > 0 and r = max{l, r}. = vr/iV’o and A)) = tt/iAq are used at the initial time step. 

The linearized backward Euler Galerkin-mixed FEM scheme (j2.10p - (|2.12l) is uncoupled and 
and (<tJ),A)() can be solved simultaneously. Moreover, for each FEM equation, one only 
needs to solve a linear system at each time step. 

Remark 2.1 A linearized Galerkin-mixed FEM in two dimensional space.' If the su¬ 
perconductor is a long cylinder in the z-direction with a finite cross section and the external 
applied field He = i7e[0,0,1]^ (i.e., He is parallel to the z-axis), the original three dimensional 
equations can he reduced to a two dimensional equation |12l [5^ 

^ - ZK(div A)fi (-V -h A)^fi -h (IV'P - l)'fi = 0 , (2.13) 

UL hi 

% 

—-Vdiv A -|- cnrl curl A H--t- IV^pA = cnrl , (2-14) 
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with boundary and initial conditions 


^—= 0, curlA = i/e, A • n = 0, on9r2x(0,r], (2-15) 

'ip{x,0) ='i/jo{x), A(x, 0) = Ao(x), in O, (2-16) 


where V’ and A = [^ 1 ,^ 2 ]'^ are scalar-valued complex funetion and vector-valued real funetion, 
respeetively. The operators div, V, curl and curl in H2.13\} - i2.16\} are defined by 


dAi dA2 ^ , 
dlV A = 1 , V = 


dx 


dy 


dtp 

dx 


dfi 

dy 


, curl A = 


dA2 dAi 


dx 


dy 


, curl 


dfi dip 
dy ’ dx 


To introduce the FEM scheme for i2.13\) - h2.10p . for a quasi-uniform triangular mesh Th with 
mesh size h = max^^^eThidiam Ox}; we denote by Vf and QP the r-th order Lagrange finite 
element subspaces ofT-L^{Ll) and respectively. We also define Q'P = Qpr\H^{Q). Let 

o o 

be the r-order Raviart-Thomas finite element subspaces o/H(div), where the case r = 0 
corresponds to the lowest order Raviart-Thomas element (3 dofs). Then, by introducing a = 
curl A, the linearized backward Euler Galerkin-mixed FEM scheme is to find ipp^ G and 
(aj^, A^) G X LJp, where ajp = TThHp on dLl, such that for n = 1,2,... , N 


{DriPPi, cu,) - m((div Ap^-^)iPPl, a;,) + ((1V + , (-V + Al-^)ojf,) 

K tx 

+ - l)iPl,iOH) = 0, (2.17) 

and 

- (curlx/i, A))) = 0, (2.18) 

(DrAl, Vh) + (div Al, div Vh) + (curl o-)), v/^) + (| A)), v/j) 

= (curli^,^v,)-^ , vv ,g tp, (2.19) 

where ip^ = TThipo and A® = tt/iAq are used at the initial time step. 


Here we focus our attention on analysis of the mixed scheme and present our main results on 
optimal error estimates in the following theorem. The proof for the problem in three dimensional 
space will be given in sections 3 and 4. The proof for the two-dimensional model can be obtained 
analogously and therefore, omitted here. Numerical simulations on a slightly different scheme 
were given in [22] . Further comparison with conventional Galerkin FEMs will be presented in the 
section O We assume that the initial-boundary value problem (|l.ll) - (|1.4p has a unique solution 
satisfying the regularity 

ip G L°°{0,T-,'H^+^),iPt e L°°{0,T-,n^+^) ,iPtt G L°^{0,T-,C^) (2.20) 

and 

r A G L°°(0, T; , A* G L°°(0, T; , Au G L°°(0, T; L^), 

I divA G L°°(0,r;H''+i),(divA)j G L°°(0,T;H''+i), (2.21) 

[ CTG L°°(0,T;H'+i),crt G L°°(0, T; H^+i), 

where I > \ depends on the regularity of the domain H. 
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Theorem 2.2 Under the assumption /12.20\) - I12.21\} . there exist two positive constants ho and 
tq such that when h < ho and r < tq, the FEM systems li2.10\) - li2.1S^) and Ii2.17\ )- li2.19\) are 
uniquely solvable and the following error estimate holds 


N 

(||« - rWh + l|AJ - A"\\l,) + T Y, IK - curlA’”|||. < C.{t^ + h^‘), ( 2 . 22 ) 

“ “ m=0 

where s = min{r + 1,/} and (7* is a positive constant independent of n, h and r. 


Remark 2.3 The above theorem shows that the convergence rate of the mixed method depends 
upon the order of the FEM spaces and also the regularity of the exact solution. In m, the 
authors proved that on a nonconvex polygon, the TDGL equations possess regularity only with 
\ ^ ^ ^ < 1? where maxjOJj denotes the largest interior angle of the polygon. More 

precisely, the two dimensional TDGL equations admit a solution with I < ^ on an L-shape 
domain. Therefore, the L^-norm convergence rate of the Galerkin-mixed method for the problem 
on an L-shape domain is 0(r +We note that the convectional Lagrange FEMs converge 
to a wrong solution due to the fact that in space A E H* with I < 1 on a nonconvex domain, 
see Example \5.2\ in section O Analysis for the TDGL equations in certain three-dimensional 
geometries was given in IMi 


In the rest part of this paper, we denote by (7 a generic positive constant and e a generic small 
positive constant, which are independent of (7*, n, h and r. We present the Gagliardo-Nirenberg 
inequality in the following lemma which will be frequently used in our proofs. 

Lemma 2.4 ( Gagliardo-Nirenberg inequality i36f ): Letu be a function defined on LI inW^ and 
d^u be any partial derivative of u of order s, then 

\\d^u\\LP < C\\d^u\\lr llwllir + C\\u\\l., 


for 0 < j < m and — < a < 1 with 

j — m, — — 


1 

p 





except 1 < r < oo and m — j — j is a non-negative integer, in which case the above estimate 
holds only for ^ < a < 1. 


3 Preliminaries 

3.1 An auxiliary elliptic problem 

We consider the elliptic boundary value problem 

curl curl u — V div u = f in H, (3-1) 

curl uxn = 0, un = 0 on dLl, (3.2) 

where u = [ui ,U 2 and f = [/i , /2 , fs]'^■ By noting the definition of the standard three 

dimensional operators curl, V and div, it is easy to verify that curl curl u — Vdivu = — Au. 
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By introducing the new variable cr = curl u, the mixed formulation of (|3.1I) with the boundary 
condition (|3.2p is given by 


cr = curl u 

in n, 


(3.3) 

curl cr — V div u = f 

in n, 


(3.4) 

crxn = 0, un = 0 

on dLl. 


(3.5) 


0 

o 



Then, the weak formulation of the above equation is to hnd (<t,u) S H(curl) x H(div), such 
that 

(o-,x) - (curlx,u) = 0, Vx:GH(curl), (3.6) 

(curl <T , v) + (div u , divv) = (f , v), V v G H(div). (3.7) 

Based on the above mixed weak formulation, the mixed FEM approximation to (I3.6p - (I3.7I) is to 

O o 

hnd (<T/j,u/i) G Qh ^ U)), such that 

(o-h,X/i) - (curlx,,,u,j) = 0, Vx/,gQ);+S (3.8) 

o 

(curl v/j) + (div u,j,divv/j) = (f,v/,), V G U)). (3.9) 


Theoretical analysis on convergence and stability of the above mixed hnite element methods 
can be found in the two seminal papers 01 [5], while our main concern in this paper is on 
the nonlinear parabolic problem of superconductivity. We summarize the main results in the 
following lemma and we refer to 01E] for details. 


Lemma 3.1 Let (<t,u) and (<t/i,u/i) be the solution of and respectively. 

Then the following error estimates hold 


- <^/i||H(curl) + l|u - U/j||H(div) 

<c( in/ ||<T - x,j||H(curi) + inf ||u - Vft||H(div)) (3.10) 

x,.6q;;+i v^eu- 


3.2 Elliptic and mixed projections 

To prove the optimal error estimates of the linearized Galerkin-mixed FEM scheme (I2.10p - 
(j2.12p . we dehne a Ritz projection operator : 'R^(n) and a mixed projection operator 

Ph ■ H(curl) X H(div) —)• x Uj) as follows (also see 00]) : for given t G [0, T], hnd Rhip G VH 

and Ph(o', A) := {P^^cr, A), P^{cr, A)) G x Uj) with P^icr, A) x n = vr/icr x n on dLl, such 
that 

- Rhi!) ,ViVh) + Mifj - Rhip,iVh) = 0, VwfeGV^ (3.11) 

where M is a positive constant which is chosen to ensure the coercivity of (13.lip and 

(curlx/j,A-P^(cr,A)) =0, Vx^gQ^+S (3.12) 

(curl(cr - P;)(cr, A)), v/j) + (div(A - P^(cr, A)) ,divv/j) = 0, V v/^ G Uj). (3.13) 


We denote the projection error functions by 

= Rhip -Ip, Pa = A) - a 


PA = Pi{o-,A) - A. 










With the regularity assumption (I2.20I) - (I2.21I) . by standard finite element theory [8] and Lemma 
EH we have the following error estimates 


< (3-14) 

I l|P<T||H{curl) < , IIPAllH(div) + \\{PA)t\\L^ < , 

and the stability result 

||-R/iV’IU°° < C, WRhipWw^’^ < C. (3.15) 


4 The proof of Theorem 12.2 


For n = 0, ..we denote 


= el = al-Pl,{a-,A^), e\ = Al - Pfi{a\ A^). 


In this section, we prove that the following inequality holds for n = 0, ..N 

r[\\eWm + Ildl22 + lldiv6^1122 } < 

m=l 


TL 

+ 116X11^2 + ^ + ||e™||^2 + ||div6^1112^ < (4-1) 


by mathematical induction. Theorem 12.21 follows immediately from the the projection error 
estimates in (j3.14p and the above inequality. 

Since 


+ II^aIIlz = IKhV’o - Rhil^^Wh + ||7r,,Ao -P^(cr°, A °)|||2 < Cih"^^, 

(BU) holds for n = 0 if we require ^ > Ci, we can assume that BH) holds for n < k — 1 for 
some k > 1. We shall find a constant C*, which is independent of n, h, r, such that (14. ip holds 
for n < k. The generic positive constant C in the rest part of this paper is independent of C*. 

From the mixed variational form (j2.6p - (j2.8p . the linearized Galerkin-mixed FEM scheme 
(j2.10ll - (j2.12l) . and the projection (I3.11I) - (I3.13I) . the error functions e^, e” and e\ satisfy 

{Drel,Uh) + 4(Ve; , Vcu,,) + M(e; 

Kj 

= -{R>Tp^,^h) + iK ((div Al-^)'il:l - (div A”"^)'0’", u^h) 

+-[Al-^rh - A"-V" , Vo;,) - - {{vrh , - (vr , A^-^uh)) 

Kj K 

+ {M^h + (1 - - \Ai-^?Wh - Mr - (1 - \r-^? -\A^-^\^)r ,u^h) 

+Rl 

5 

= Y,J?{^h) + Rl. Viv^evl n = l,2,...,N (4.2) 
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and 


h > n = 1,2,... 


where 


and 


{e^^Xh) - {curlxh,e\) = 0, Vxfe€Q''+^ 

{Dre\,Wh) + (dive^ ,divv/j) + (curie" ,v/j) 

- {{r-^rvr-^ -r-^v{r-^r,^r,,)] 

8 

= + Vv;,gUL n = l,2,...,iV, 

i=Q 


= {Drr - + in (div (A" - A^-^r ,^h) 

+ -{Vr , (A"-i - A")a;0 - -((A"-i - A^)^ , Vo;;,) 

K K 

+ ((|A"-Y - |A"p + IV’"-^ - m^)y",ujh) 


R\ = {DrA^ - - {\r\^A^ - IV^^-^A",V,) 


(4.3) 


(4.4) 


dt 


--iurr^r-r^irr) - ar-r^r-^ -r-^^ir-r)Ah) 

define the truncation errors. 

We take uih = e"^ and {xh ,'^h) = (e^ >6^) in (I4.2p - (j4.4h . respectively. It is easy to see that, 
with the regularity assumptions ()2.20p - ()2.2ip . 


|i?;| + \Rl\ < e||Ve;||i. + C7||e;||i. + C||el||i. + 


(4.5) 


We now estimate for i = 1, ..., 5 in (14.21) and J"(e^) for f = 6, 7,8 in (|4.4p term by 

term. By noting the projection error estimates (I3.14p . we have 

|jr(e;)l<C||e;f+ Ch2% 

\J^{e\)\<C\\e\f + Ch^^ 


and 


^(-^2"(e;)) = 7ec(m((divA"-i)e; + (divA"-i)i?;,V'”-(divA"-i)V’",e;)) 
= Re (iK((div A"-^)i?,,V’” - (div A"-i)V^" , e^)) 

< C|((div(er' + /^r'))^/^V’",e;)|+C|((divA"-^)p;,e;)| 

< C{\\dAeXy\L2 + ||divpl-i||^3)p,^"|U^||e;||i. 
+C||divA"-i||i^||p;-i||^.||e;|U. 

< e||diveX-i||i2+e-^C(||e;||i.+h2^), 

where 77e(Ji) denotes the real part of R and also, we have noted the fact that 

77£((iK(div A""^)e^,e();)) =0. 
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For the term Jg , we see that 

|7^(J3"(e;))| < C|Xm(((er' + pVWh + A"-'(e; + , Ve;))| 

< C|Xm(((er' + Ve;))| + C\Tm[{{e^-^ + , Ve^)) | 

+ C|Xm((A-i(e; + 4),Ve;))| 

< C|Xm((el-ie;^,Ve;))l+ C||pX-ii2||e;|U^||Ve;|U2+e||Ve;||i. 

+ e-'C(||e;||i. + ||el-'||i.+h 2 ^) 

< C|Xm((er'e;,Ve;))| +Ch-5||e;||^, +e||e;||^i 

+ e-'C(||e;||i. + ||er'lli 2 + /.^^) 

< C|Xm((er'e; , Ve;))| + e||e;||^i + e-^C(||e;||i. + Ijer'lli^ + 

where we have used the projection error estimate (13.141) and an inverse inequality and required 
< e. 

Similarly, J 4 is bounded by 

\MJi{eV)\ < C|Im((V(e; + 4), A;-‘e;)) I + C\Im((Vr , + Pr‘)4))l 

< C|Im{(V(eJ + pD , (ey■ + | + C|Im((V(4 + pj), A”-'^)) | 

+f|| 65 ll^i + f ^C'(||e ^||^2 + \\e\ ^ 11^2 + /i^®) ■ 

By requiring Ch ^~2 < e and using inverse inequalities and integration by parts, 

|Xm(Vp;,(er'+Pr')«;)) I < l|Vp;||i3||er'+Pr'llL3||e;ilLa 

< C(h* + ||er'llL^)l|e;iki 

< 4e^4m+C{\\el-YL^+h^n, 

|Xm((Ve:^,pX-ie;))| < ||Ve;^||^2 IIp^'IIl^ lle^lk- < ||e;||l,i < e||e;|||,i , 


|Xm((V(e;^ + p;),A’^-ie:^))| = |Xm((e; + , (div A^'^je;^)) | 

+ |Xm((e; + p;,A"-i.Ve;))| 

< e||e;||^,+e-^C||e;||i 2 +e-W^ 

Then it follows that 

|7ec(jr(e;))l < C|Xm((Ve;,er'e;))| + 6 ||e;||^, + 6 -^C(||e;||i 2 + ||el-ii 2 + 
For the term J 5 , we have the bound 

-ReiJ^ie^)) < -7^e((K■l|2 + |A^ll 2 )V.;:-(|V^"-Y + |A"-l| 2 )^^e;) 

+(M + l)||e;||i 2 + (M + l)(p;,e;). 

Since 


-7J<;{(|Ar>|2«-|A'‘-‘|4”.4)) 

= -(|Ar'l"4, ej) - 7fe((|Ari|2flrf" - |A”-1|4” , 4 )) 

< -(lArTej, ej) + |(|A”-i|Vj ,4)1 + |((|a;-‘I" - | A”->4/;rf” , 4)1 

< -(I Ar T4. 4) + l((Ar' + A”->). (el-' + pl-')Rer ,4)1 

+ e||e^||^i+e ^Ch?'^ 

<6||e;||l,.+6-iC(||er'lli. + /^^*) 
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and 


the real part of can be bounded by 

neiJ^e-^)) < e||e;||^,+6-iC(||e;-i||i. + ||er'lli2+/i"^). 

For Jy , we can see that 

j^ei) < c|xm(((v^r')*v^r 1 - ir-^yvr-^, cd) \ 
+c|xm((v^r'v(V'r')* - r-^vir-^r,ei))\. 


(4.6) 


Moreover, the first term in the right hand side of (14.61) is bounded by 

ii™((«'r‘)'vvr‘ - .ei))i 

= |im(((«-‘rv(e;->+p;-') + (e;-'+p;-'rv^”-‘,ei))i 
<C|Im(((V.r‘rv(e;->+p;-‘),el))t+t||e;->||J,,+e-‘C(||el||J. + ft"*) 

<C|Im(((e;-‘rv(e;-‘+p;-‘),el))| + C|Im{((Rrf”-'rV(e;-‘+p;-‘), 4 ))| 

+ 4e"^-VH.+‘-'C{\\el\\l, + h'^‘) 

+ C|Im(((ii,.r-‘rVp;->.el))| +e||e;-‘||^. +e->C(||el||i. + ||e;->||i. +/,."*). 

< C|Im(((e;-‘)*Ve;-> ,el))| +C/»-i||e;->||„.||ei||y 

+ e||div6X11^2 + e||e^ ^||^i+e ^^(116X11^2 + 116)^; ^||^2 + /i^*). 

< C|Xm(((e;-i)*Ve;-SeX))| +e||diveX||i2 +6||e;-i||^i 
+ e-^C{\\el\\l, + \\e-fX^ + h^y, 

where we have used an inverse inequality and required that Ch^~^ < 1. Similarly, the second 
term in the right hand side of (|4.6p is bounded by 

c\imM-^v{rH-^r - r-^v{r-^r , eX)) i 
< C|Xm((e;-^V(e;-i)* ,eX))| + 6 ||diveX||i 2 + 6 ||e;-i||^, 

+ e-^C(||eXlli2 + ||6;-i||i2+/.2*), 

It follows that 


J^el) < C\Im{{{e; 


n— 1 \* 


Ve; 


n—1 


eX))l + f||diveX||i2 + e||e” ^||^i 


+e-^C{\\el\\i, + \\e;-X,+h^y. 

Finally, for the term J^, we can derive that 

JSiel) < -(K■TA5(-|V^”-l|2A^eX) + C||eX||i 2 +C'/. 2 ^ 


= -(IV’r'p6X,6X) - (IV’rY^h^K, A") - \r-^\^Al,el) + C\\el\\l,+Ch^^ 
< -dV^r'PeX ,6X) + C|((V^ri + V^"-^)(e;-^ + p;-^)p2(^-, A") ,eX)| 
+C\\el\\h+Ch^^ 


< elle; 


n-l ||2 


■j^,+e-^C{\\el\\l, + \\el-y\l,+fi^n. 


(4.7) 
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With (14.51) and above estimates for Jj, for z = 1,..., 8, adding (14.3p . (14.41) and the real part 
of (j4.2l) together, we get 

'Re{{DTe^,e^)) + {Dre\,e\) + ■^||Ve^||^2 + ||diveX||i2 + ||e”||^2 

rv 

< C|Xm(((e;-^)*Ve;-Sel))| + C|Xm((er, Ve^^)) | 

+ 6 (||divel||i 2 + lldiver'lli^ + \K\\h^ + 

+e-^C{\\e-4l, + \\e;-X^ + \\el\\l, + , (4.8) 

for n = 1, ..., A;. 

In order to estimate |2m(((e^“^)*Ve^“^, e^)) | and |2m((e^“^e();, Ve();))| in (14.Sp . we use the 
induction assumption ()4.ip for n < A: — 1 and inverse inequalities. 

For r < /i®, we have 


|I™(((er‘)'^«r‘>«A))l < C||Ve;-‘||i.||e;-‘||i.||el||t3 

< C||e;-‘||?,,/r4||el||i.. 

To estimate ||6^11^2, from (14.41) and estimates (14.7p . we have that 


l|ellli 2 +T(||e-||i 2 + ||div 6^11^2) 

+ r{e\\e;-YH^+e-^Ci\\el\\l, + We^-^L^ + + h^^)) 

< (-+ C'r)||eX|||2 + ||e^ ^11^2 + T^CIIeJJ; ^W^ih ^ 

<l\\el\\h+Cah^\ 


where Ct < ^. Then substituting the last estimate into (14.9p gives 


where we require Cy/Clh ^~"2 < e. Similarly, 


< 4eT^fm 


C|Xm((er'e;,Ve;))| < 

< 

< 

< 

< 

< 


^^l|er'llL3||Ve;|U2||e;||i6 

t^l|er'llL3||Ve;^||i2 

Ch-"^\\eV\\L4^e^\\h 

ch-^y^e^^ 7 h^||ve;iii 2 

Cv^h*-i||Ve;||i 2 


where h satisfies Cy/C^h^ 


2 


< e. 


(4.9) 


(4.10) 


(4.11) 


13 



For /i* < T, by (14.11) . we have that for n <k — 1 

Hernia + 116X11^2 < , 

IIc^IIhi + II^ctIIlz + ||div eW\\2 < T < C:tT . 

Then, we have 


C|Xm(((e;-i)*Ve;-\el))| < 

< 

< 


C’l|elllL3||e;-lHi||e;-'ka 
Cv^(||e^||i 2 + ||divel||i2)||e;-i||Hi 
e(||e^||i. + ||divel||i. + ||e;-i||2^0 


where r satisfies CC^t < e and we used the following discrete embedding inequality (the proof 
is given in appendix) 


l|elllL 3 < C(||e"||i. + ||divel||i.) 

By noting (jd.lip and (I4.12p . we have 

C|Im((el->e5,V4))l < C||4-> b, ||4||?,. 


(4.12) 


A 

,n—1 

cr 


< C{\\e 

/C 

< 4e^4H^ 


< Cv^r||e ""2 


IIl 2 + lldive^ 

ipwm 


-\\l^)\K\\ 


m 


where we require that r satisfies Cy/C^T < e. 

Therefore, from (14.Sp . for both cases t < and < r, we can derive that 

7^((Zl.e;,e;)) + {Drel,el) + ^||Ve;||i. + M||e;||i. + ||divel||i. + \K\\l, 
< 6(||e^||i. + lldivellli, + llejlll,^ + ||e;-^||l,0 
+e-iC(||e;||i. + \\e;-X^ + \\el\\l, + ||erii2 + . 

By choosing a small e and summing up the last inequality from 1 to n, we arrive at 


n 

114111= + Il4lli= + T {|i4Pff, + WdiveXWh + Il4lli=) 

m=l 

n 

<CtY. (Il4lli= + l|6Alli=) + C(t2 + /,2') . 

m=l 

With the help of the Gronwall’s inequality, we can deduce that 


Il4llh + Il4llh + E ^(ll4llb + licili= + l|3iveSiii.) < C(4 + h^-). 

m=l 


Thus (|4.1I) holds for n = /c, if we take C* > 2C. The induction is closed and the proof of 
Theorem 12.21 is complete. I 
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5 Numerical results 


In this section, we provide several numerical examples in both two and three dimensional spaces 
to confirm our theoretical analysis and show the efficiency of the linearized Galerkin-mixed FEM 
scheme. The computations are carried out with the free software FEniCS [33]. 


5.1 Two dimensional numerical experiments 


Example 5.1 In this example, we consider the following two dimensional artificial problem 


^ - iK(div A)ip + (-V + A) V + = S' , (5-1) 

Ot K 

—- Vdiv A + curl curl A H - — 'ijj’Vijj*) + I'iApA = curl + i , (5.2) 

(jt 2 /^ 

for t G (0,T], X G II, with boundary and initial conditions 

dib 

= 0, curl A = Hf., A • n = 0, on dll, 

dn 

V’(x, 0) = V’o(3;)) A(x, 0) = Ao(x), in n. 


where II = (0,1) x (0,1) and k = 1. The functions f, g, i/jq and Aq are chosen correspondingly 
to the exact solution 


V’= exp (—l)(cos(7rx) + i cos(7ry)), A = 


exp(y — t) sin(7rx) 
exp(x — t) sin(7ry) 


with 


He = exp(x — t) sin(7ry) — exp(y — t) sin(7rx). 


We set T = 1.0 in this example. 



X 

Figure 1: A uniform mesh on the unit square with M = 10 (Example [5T]). 

We use a uniform triangular partition with M + 1 vertices in each direction, see Figured) for 
M = 10, where h = We solve the system (I5.1|) - (I5.2|1 by the proposed linearized backward 
Euler Galerkin-mixed finite element methods (I2.17l) - (|2.19l) with (V^, x U))) for r = 0, 1, 
2 . As c 7 = curl A is a real scalar function in the two dimensional case, correspondingly its 
approximation space is the Lagrange element space consisting of all piecewise polynomials 
up to (r + l)-th order. As the exact solution is smooth,to confirm our L^-norm error estimates. 


15 



























Table 1: L^-norm errors of A and a on the unit square fExaiuple 15.11) .. 


r = 

1 

M 


Ajy - A(-,f 7 v) L 2 

hh - (^{■,iN)\\L^ 

M = 64 


3.1216e-02 

1.0296e-02 

2.0804e-03 

M = 128 


1.6284e-02 

5.1766e-03 

1.0637e-03 

M = 256 


8.3156e-03 

2.5959e-03 

5.3770e-04 

order 

0.95 

0.99 

0.98 

(vtQkXu;.) T = 

1 

TP 


W^h ~ A(-,tAf)||L2 

hh - (^{■,iN)\\L^ 

M = 16 


9.9391e-03 

8.1354e-04 

6.3113e-04 

M = 32 


2.4956e-03 

2.0431e-04 

1.5807e-04 

M = 64: 


6.2454e-04 

5.1132e-05 

3.9533e-05 

order 

2.00 

2.00 

2.00 

(Vj.QjxUj) T = 

1 

TP 


- A(-,f7v) l 2 

\Wh - 

M = 8 


4.1680e-03 

4.5426e-04 

2.7586e-04 

M = 16 


5.2687e-04 

5.7080e-05 

3.4243e-05 

M = 82 


6.6130e-05 

7.1400e-06 

4.2634e-06 

order 

2.99 

3.00 

3.01 


we set r = ( 7 ^)*'"*"^ in our simulation and we present in Table [U the L^-norm errors of "0, A, a. 
We can see clearly that the L^-norm errors are proportional to r = 0, 1, 2, which confirm 

the optimal convergence of the scheme in two dimensional space. 

To show the unconditional convergence of the methods, we solve (I5.1I) - (I5.2I) by the scheme 
(l2l71) - (l2Tm with (V^, X Uj,) and three different time steps r = 0.1, 0.01, 0.001. For each 
fixed r, we take M = 8 , 16, 32, 64, 128. The errors of ?/>, A and a are presented in Figure 
[2j From Figure O we can see that for each fixed r, when the mesh is refined gradually, each 
error converges to a small constant of 0(t), which shows clearly that the proposed scheme 
(I2.17p - (I2.19I) is unconditionally stable and no time step restriction is needed. 





Figure 2: errors of A and a with x Uj^) iFxample 15.11) . 
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Example 5.2 In this example, we investigate the TDGL equations (I5.1l) - (l5.2p on an L-shape 
domain with a nonsmooth solution by the Galerkin-mixed method. This example was studied in 
m by a projected FEM based on Hodge decomposition. Here, we use the same exact solution as 
in m, i e., K = 10, and the functions f, g, Hg = curl A, i/jq and Aq are chosen correspondingly 
to the following 'ip and A 


i’ 


t^<h(r)r3 cos 



3 + 3 ^ pQg 

(t)- 

^|f^<I>(r)r“3 + t2<I)'(r)r3^ sin (|) 


where {r,9) denotes the two-dimensional polar coordinates. <h(s) in the above expressions is a 
cut-off function defined by 


$(s) 


0.1 ifs<0 .1, 

T(s) if0.1<s<0.4, 
0 if s > 0.4, 


where T(s) is a seventh order polynomial satisfying 

' T(O.l) = 0.1, T(0.4) = 0, 

< T(^)(0.1) = T(2)(0.1) = T®(0.1) = 0, 

T(^)(0.4) = T(2)(0.4) = T(3)(0.4) = 0. 

It is noted that in spatial space, the above exact solution only has the regularity 
V'e "Ht'S A G Hi“^ (T G divAG^^t“^ for e > 0. 



Figure 3: A uniform mesh on the L-shape domain with M = 4 (Example 15.21) . 

Now, we use the Galerkin-mixed method (j2.17p - (j2.19jl to solve the above problem on uniform 
meshes as shown in Figure [3l where h = ^. We set the the terminal time T =1.0 and do 
computation with {'iph,Ah,ah) G (V)^,Qjj,U°) and the time step t = ^. The L^ errors are 
shown in Table [5J From Table [21 we can see that the convergence rates of the Galerkin-mixed 
method for the three components pj, A and a, are in the order of and 

respectively, which indicates that our method is suitable for problems on nonconvex domains. 

For comparison, we also solve this artificial problem by a Lagrange type FEM under the same 
data settings. The linearized backward Euler Lagrange FEM is to find V’/J G Vp and A^J G V^, 
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Table 2: norm errors of 'i/’, A and a by (V^, x U^) on the L-shape domain fExamole 15.21) . 


r = ± 

IIV'^ - V'(-,%)IIl2 

W-^h ~ A(-,t7v) Ll 2 

hh - cri-,tN)\\L^ 

M = 32 

6.5548e-04 

1.7728e-02 

1.2104e-01 

M = 64 

2.6140e-04 

9.2365e-03 

3.3953e-02 

M = 128 

1.2130e-04 

5.0878e-03 

8.9061e-03 

M = 256 

5.9174e-05 

2.9632e-03 

2.2696e-03 

order 

1.04 

0.86 

1.91 


such that for n = 1, 2 ,... , 

{Drrn.^h) - m((div AI-^Wh.u:h)+{C-V + AI-^Wh , (^V + Ari)a;,) 

-((lAr^P + = 0 Vcu, G V^ (5.3) 

and 

{DrAl,Wh) + (divA^divvO + (curl A)), curl V;,) + ^A)), v;,) 

= (if-,curlv,)-^((V^r')*VV^ri-V^rV(^r')*,v,) Vv.gVI, (5.4) 

where denotes the linear Lagrange FE subspace of {A | A G H^, A • = O}- V'h = tt/iV'o 

and A^ = tt/jAo are used at the initial time step. The errors obtained by the conventional 
Lagrange FEM (|5.3l) - (|5.4p with M = 32, 64, 128, 256 are presented in Table El We see clearly 
that the numerical solution of the Lagrange FEM does not converge. We refer to [221 [32] for 
the numerical observation of the spurious convergence phenomenon by conventional Lagrange 
FEMs in the vortex motion simulations on an L-shape domain. 

Table 3: L^ norm errors of A and a on the L-shape domain by the conventional Lagrange 
FEM (|5.3|) - (|5.4I) lExamole 15.21) . 


ivlAl) 

r = A 

\\'>Ph - V'(-,i7v)||L2 

A^ - A(-,fAr) ^2 

- 0-(-Tjv)IIl2 

M = 

32 

3.8735e-03 

9.6105e-02 

5.1970e-01 

M = 

64 

3.8616e-03 

9.8059e-02 

4.7659e-01 

M = 

128 

3.8466e-03 

1.0298e-01 

4.6553e-01 

M = 

256 

3.8395e-03 

1.0539e-01 

4.6204e-01 


5.2 Three dimensional numerical experiments 

Numerical results for a three-dimensional artificial example was reported in [22]. For complete¬ 
ness, we include some similar results in this section to confirm the convergence rate of the 
Galerkin-mixed method for the problems in three-dimensional space. 
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Example 5.3 We consider the following three dimensional artificial problem 


^ - iK(div A)V’ + (^V + + (IV’P = g, 

^ yV 7 

— -Vdiv A + curl curl A H-+ IV’pA = curlHg + f , 

(yt> 

for t E (0,T], X E il, with 



V^(x,0) = 


curl A X n = Hg x n, A • n = 0, 

'^o{x), A(x,0) = Ao(x), 


on (9fl, 
in fl, 


(5.5) 

(5.6) 


where we set n = (0,1) x (0,1) x (0,1) and k = 1. The functions f, g, ipo and Aq are chosen 
correspondingly to the exact solution 


V' = exp (t) (cos(7rx) cos(7rz) + i cos(7ry) cos(7rz)) , A = 


exp(t) sin(27rx) sin(27ry) 
exp(t) sin(27r?/) sin(27rz) 
exp(t) sin(27r2;) 


with 


He 


—27r exp(t) sin(27ry) cos(27r2;) 
0 

—27r exp(t) sin(27r3:) cos(27r?/) 


We also set the terminate time T = 1.0 in this example. 

A uniform tetrahedral mesh with M + 1 vertices in each direction of the cube is used in our 
computation, where h = ^. We solve the system by the Galerkin-mixed FEM scheme (I2.10p - 
(I2T^ with (V^, X U))) for r = 0, 1, respectively. Here cr = curl A is a three dimensional 
vector function and correspondingly is the (r+l)-th order first type NMelec element space. 
To confirm our error estimates in the norm, we choose r = for (V^, x u;;). We 

present in Table|l]the L^-norm errors of ip, A, a. We can see clearly that the norm errors are 
in the order of r = 0, 1, which confirm our error analysis in three dimensional space. 


Table 4: norm errors of ip, A and cr on the unit cube (Example [5]3|). 



Ql 

X UO) T = 

1 

M 

IIV’^ - V’(-Tiv)llL2 

W-^h - A(-,tAr)||L2 

||cr^ - cr(-,t7v)||L2 



00 

II 


2.0951e-01 

7.0869e-01 

3.9600e+00 



M = 16 


8.3881e-02 

3.3623e-01 

1.9502e+00 



M = 32 


3.7858e-02 

1.6561e-01 

9.7032e-01 

order 

1.234 

1.049 

1.015 


Ql 

X U),) r = 

1 

JF 

IIV’^ - V’(-Tiv)llL2 

W^h - A(-,tAr)||L2 

||cr^ - cr(-,t7v)||L2 



II 


4.2803e-01 

3.3717e-01 

1.6782e+00 



00 

II 


1.0461e-01 

8.9300e-02 

4.4795e-01 



M = 16 


2.5907e-02 

2.2832e-02 

1.1670e-01 

order 

2.023 

1.942 

1.923 
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6 Conclusions 


We analyzed a linearized backward Euler Galerkin-mixed FEM for the time-dependent Ginzburg- 
Landau equations in both two and three dimensional spaces. We have established uncondition¬ 
ally optimal error estimates for the three dimensional model, while analysis presented in this 
paper can be extended to many other cases, such as the two dimensional problem and Grank- 
Nicolson scheme. The method can solve for the induced magnetic field curl A (curl A in two 
dimensional space) accurately without corner singularities and the method is particularly suit¬ 
able for nonconvex domains and domains with complex geometries. Numerical experiments 
presented in this paper show the efficiency of the method and confirm our theoretical analysis. 
Large scale parallel computations and adaptive local refinement or moving mesh strategy will 
be conducted in our future work. 
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Appendix 

o o 1 

With all the notations in section [21 for any given e\ G if there exists a e” G such that 

(e^,Xh) - (curlx/„eX) = 0, VG (6.1) 

then the following discrete Sobolev embedding inequality holds 


IcaIIls < C* (116^1^2 + ||diveA||L2) 


( 6 . 2 ) 


Proof: 

From the Hodge decomposition HE], the following exact sequence holds: 


H(curl) - 

H(div) ^ 



i 


0G- - 

curl -^-rr div 

—^ ^ 

o 

cr 

where Lq = {u G {v, 1) = 0} and S)) 

= S)) n Lq with 



SI :={v e L\n)\ynK eTh, v\nii is a polynomial of degree r}. 

Then, we have 


UD = curl Q 


r+l 


grad^jS;;, 


(6.3) 


where 0 denotes the direct sum and the linear operator grad;^ : S)) —)• U)) is dehned as: for any 

o o 

given Sh G S)), hnd grad^^s/i G U)) such that 


{gradf^Sh,Vh) = -{sh,divvh), Vvh G U)). 


Therefore, for any given G U)), there exist 6^ G and Sh G S)) such that 


To prove 


Ca = curl^fe 0 grad/jSfe 
we only need to show 

11 grad;, s/j 11^3 0 11 curl 6 »/, 11^3 < C (||e” 11^,2 0 ||div 6 ^ 11 ^ 2 ) . 


(6.4) 


(6.5) 


We first estimate grad;,s/i. Indeed, grad;,s/, can be viewed as the mixed FEM solution to 
the following Poisson equation with pure Neumann boundary condition 


—A u = div Ca in ^ 


Vu • n = 0 


on 


on 


( 6 . 6 ) 


It should be noted that the regularity for (|6.6p depends upon the the domain H, see usmu. 
For a general polyhedron in three dimensional space, we have the following shift theorem 


\\u\\jj3/2+c, < C\\div eX\\H-i/ 2 +o , (6.7) 

where a > 0 depends only upon the geometry of the polyhedron. The classical mixed FEM for 
solving (j6.6p is to find {qh,Uh) S (U)), S))) such that 

I {qh,Vh) + {uh,divvh) = 0 Vu;,gU)), 

1 -(div qh, iJ.h) = (div e%, fit) V /r/* G S)). 
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From ()6.4I) and (16.5p . one can verify that grad/jS/i = qh. By using an inverse inequality and the 
classical error estimates of mixed methods for elliptic equation mi, we can deduce that 

11 grad;, Sft 11^3 < ||grad;,s;i - 7r;,Vu||i3 + ||7r;,Vu||i3 

< Ch ~2 ||grad;,s/, - 7r/,Vu||2,2 + ||u|| jj3/2+c. 

< C/i"2(||grad;,Sft - Vu\\l2 + ||7r;,Vu - Vu\\l2) + CUdiv6X1^2 

< C7/i-U^+"||Vu||j^i/ 2+. +C||divel||i 2 

< C||div6^11^2 , (6.8) 

where -Kh is the projection operator defined in (12.91) . 

Now we turn to estimate curl0/,. From the discrete Hodge decomposition (16.Sp . we have 

o ... 

(curl Gh, curl Xh) = (el, curl Xh) = {e^,Xh) for any Xh G Ql ■ 

A key observation is that 0/, can be viewed as the mixed finite element solution to the vector 
Poisson equation 


curl curl u — V div u = e” in H 


M X n = 0, divu = 0, 


on dQ 


(6.9) 


o o .. o 

which is to find {4>h,Uh) G V))) x , where Vj^ denote the r-th Lagrange element space with 
zero trace, such that 


-{4>h,Ch) + {uhVCh) = 0 , 


for any (h G V)[. 


\ C^(l>h,Xh) + {curl Uh, curl Xh) = ie'^,Xh), for any e Q);+^ 

Since G^ = u^, from the standard error estimate mis], we have 

IldivU - 4>h\\L2 + \\U - 0h||H(curl) 

<c( inf lldivu - C/i||l 2 + inf ||u - x,,||H(curi)) 

Thus, with the projection operator tt/, in p2.9p . we can derive the following estimates 


( 6 . 10 ) 


( 6 . 11 ) 


|curl0/i||£3 < ||curl0/i — curl TT/iU 112,3 + II curl 7r;,M||j;^3 

< Ch “2 ||curl 0 /, — curlu||2,2 + ||curl 7 r/,u — curlu||2,2 + ||e”||2,2 

< C/i“2 llcurlvT/jU — curlu||2,2 + ||e” ||2,2 

< {\\curlu\\^i/ 2 +a + ||div u|| 22 i/ 2 +c:.) + ||e”||i 2 

< C||e^|| 2 , 2 , 


where we have used an inverse inequality and noted that in p 6 . 10 l) 

||curlu|| 22 i/ 2 +c, + IIdivu|| 221 / 2+0 < C'||e”|| 2,2 


( 6 . 12 ) 


(6.13) 


for a certain a > 0 , see mm- 

The embedding inequality (| 6 . 2 p is proved by combining 


and (| 6 . 12 p . 
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Remark 6.1 The discrete embedding inequality h6.0j} is enough in the error analysis of the 
Galerkin-mixed FEM for the TDGL equations. However, if the polyhedron Tl is convex we can 
derive a stronger discrete embedding inequality 


||eXllL6 < C die" 11^2 + ||diveX||L2) . (6.14) 

On a convex domain, we have a > ^ in for the shift theorem. And then, we have 

||grad;,s,j||i6 < ||grad;,s/i - TThVu\\L6 + llvr/iVull^e 

< Ch~^\\gradf,^sh - TTh'^u\\L2 + \\u\\h2 

< C||div 112,2 . (6.15) 

Moreover, for a convex domain Q, the embedding inequality i6.13\) can be re-written by 

||curlu||j2i + ||divu||22i < C'||e"||2,2, (6.16) 


which implies 


I|curl0,,||2,6 < C||e" 112,2. (6.17) 

Remark 6.2 The discrete Hodge decomposition < I6'.5'|) is only valid for simply-connected domain. 
If the domain H is multi-connected, we have 

U); = curl © grad^jS); © Har)), (6.18) 

o o .j 

where Har)) := {v^ e U)) | divvh = 0, {v^ ,curlx/i) = 0,'^Xh G i stands for the harmonic 

o 

part o/U)). In this case, 

e\ = curl Gh © grad^^s/, © Vh , 

where Vh G Har)). we can prove the following inequality via a similar analysis 

< c die" 112,2 + ||diveX||L 2 ) • 
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